Predicting the self-assembly of a model colloidal crystal 
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We investigate the self-assembly (crystallisation) of particles with hard cores and isotropic, square- 
well interactions, using a Monte Carlo scheme to simulate overdamped Langevin dynamics. We 
measure correlation and response functions during the early stages of assembly, and we analyse the 
results using fluctuation-dissipation theorems, aiming to predict which systems will self-assemble 
successfully and which will get stuck in disordered states. The early-time correlation and response 
measurements are made before significant crystallisation has taken place, indicating that dynamical 
measurements are valuable in measuring a system's propensity for kinetic trapping. 



I. INTRODUCTION 

The self-assembly of individual components into or- 
dered structures has been studied in a variety of contexts 
including biomaterials research [Ij [5] and nanoengineer- 
ing [3l-[9]. Significant progress has been made in exper- 
imentally synthesising different building blocks designed 
to assemble into specific target products [TUl - fT^ . As- 
sembling systems, however, often get stuck in metastable 
disordered states before reaching the equilibrium ordered 
ones. Kinetic traps are absent from classical theories of 
phase change [13j but they pose real problems in exper- 
iments such as those on protein and colloidal crystalli- 
sation [I4l - [l6] . In order to address the kinetics of self- 
assembly, out-of-equilibrium theories need to be devel- 
oped. 

Crystallisation is a vifell-studied process, whereby par- 
ticles in a fluid state undergo a phase transition and de- 
velop long-ranged crystalline order. This spontaneous 
ordering is an example of self-assembly. Controlling the 
crystallisation of proteins and colloids would have appli- 
cations in biology (for example, in determining the struc- 
ture of biomolecules [17] ) as well as in photonics (for ex- 
ample, in controlling the propagation of light [THIIIS]). In 
both colloidal and certain biomolecular systems the con- 
stituent particles interact via an effective short-ranged 
attractive potential with repulsive 'hard' cores [20lI26j . 

In this paper, we study a dilute suspension of attrac- 
tive colloidal particles that separate into crystal and fluid 
phases, interpreting crystallisation as a simple example of 
self-assembly, as in Ref. [23 • A schematic phase diagram 
is shown in Fig. [l]^a) , with points indicating seven char- 
acteristic bond strengths that we have considered. The 
systems we investigate begin in a far-from equilibrium 
fluid state, and evolve towards the stable crystal: we are 
interested in the dynamics of this process. The questions 
we address are the following: what happens during the 
process of assembly? How does a system that will even- 
tually assemble into a crystal differ from one that will get 
kinetically trapped? What are the first signs of frustra- 
tion and when do they appear? How does the strength 
of attractive interactions influence the dynamics? What 
are the relevant physical quantities one should measure 
in order to understand the dynamics, and at what times 
should we measure them? 
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FIG. 1. (a) Schematic phase diagram for a square-well poten- 
tial, following Ref. [28]. The legend shows the bond strengths 
U that will be considered in this work, and we show estimates 
of the locations of these state points within the schematic 
phase diagram, (b) We plot a measure of the crystallinity of 
the system, for various bond strengths, at a time t = 10® MC 
steps. The crystallinity is measured using a common neigh- 
bour analysis [29], as described in the main text. The colours 
of the symbols indicate the transition from high temperatures 
(red, weak bonds) towards low temperatures (violet, strong 
bonds). 



Self-assembly processes require a balance between a 
net drive to assembly and kinetic accessibility. The 
first requirement implies that the thermodynamic equi- 
librium state is an assembled one, which presupposes 
sufficiently strong attraction between individual compo- 
nents. The second requirement is kinetic and ensures 
that the thermodynamic state is accessible within rea- 
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sonable timescales. The two requirements are competing, 
with the thermodynamic one favouring low temperatures 
and strong bonds and the kinetic one high temperatures 
and weaker bonds. Fig. [l|b) shows the crystalhnity of 
the system as a function of the bond strength; its non- 
monotonic form demonstrates the competition between 
kinetic and thermodynamic effects, with a clear maxi- 
mum at intermediate bond strengths. 

In the past, the two competing requirements for ef- 
fective self-assembly have been characterised through 
fluctuation-dissipation ratios (FDRs), which measure the 
extent to which equilibrium fluctuation-dissipation theo- 
rems (FDTs) continue to hold in out-of-equilibrium set- 
tings. In particular, it was shown [30] that correlation 
and response functions can provide a measure of the bal- 
ance between the net drive to assembly and kinetic ac- 
cessibility, and therefore indicate which systems will as- 
semble into ordered structures and which ones will get 
kinetically trapped. Jack, Hagan and Chandler jSOj in- 
vestigated and compared two models: viral capsids and 
sticky disks. The purpose of this paper is to extend these 
ideas into a model colloidal system in order (i) to study 
whether the fate of this system can be predicted from 
the early stages of assembly based on FDR measure- 
ments, thus providing insight and subsequently guide- 
lines for colloidal and protein crystallisation experiments 
and simulations; and (ii) to identify generic features of 
self-assembly that are not system-dependent. 

Details of the computational model are given in sec- 
tion [nj In Section |III| we show how the system evolves 
in time for different interaction strengths, identifying 
regimes of slow nucleation, rapid assembly and kinetic 
trapping - features that are not predicted by the equilib- 
rium phase diagram. We then introduce correlation and 
response functions in Section [TV] explaining how they can 
be used to predict assembly quality. We include a dis- 
cussion of the robustness of our measurements: how they 
depend on measurement times, interaction range and vol- 
ume fraction. Finally, Section [V] gives a discussion of the 
results and poses questions for future investigation. 



II. MODEL 

We study a system of N spherical particles in a cu- 
bic box of volume V, with periodic boundary conditions. 
The particles have hard cores and isotropic, square-well 
interactions. Their hard-core diameter is a, the depth 
of the potential is U and the range of interaction is ^a. 
Thus, the energy of the system is 



for example in Liu et al. [5S], the short range of the at- 
tractions reduces the stability of the liquid phase, and 
the vapour-liquid binodal is metastable, lying within the 
vapour-solid coexistence curve. The phase diagram is 
sketched in Fig. [ija), where we have labelled the state 
points that we focus on in this article. The corresponding 
bond strengths are shown in the legend; for simplicity we 
have defined u = U /T which will be used to quote the 
bond strengths in the rest of the paper (we take Boltz- 
mann's constant — \ throughout). The colours of 
the symbols indicate the transition from high tempera- 
tures (red circles) to low temperatures (violet stars), in 
an order that follows the spectrum of visible light. 

We have carried out dynamical simulations of this 
model. Beginning with equilibrated systems of hard 
spheres, we use a Monte Carlo (MC) scheme as an ap- 
proximate method to simulate the overdamped Langevin 
equation 
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where E is the total energy of the system, = 
(gf"'Z^'Z^) ^ usual[3ni, and the components of the 
vectors rji are independent white noises. 

On each move of our Monte Carlo scheme, we pick 
a random particle and propose a random displacement 
from a cube of side 2ao, centred at the origin. (Thus, 
the maximum displacement in each of the x, y and z 
directions is ag.) We accept the move with probability 
min(l, e"^'^/''^'^) where A_B is the energy difference be- 
tween the states before and after the proposed move. A 
Monte Carlo step (or sweep) consists of N Monte Carlo 
moves, and we associate it with a time increment tq. In 
the limit of ao,To — >■ while holding D = o^/Qtq con- 
stant, this method provides dynamical trajectories in ac- 
cordance with the Langevin equation ([2| above [34j- In 
the Langevin description, the natural unit of time is the 
Brownian time tb = cP' jD = {a j a^^'^T^). The behaviour 
of the model depends on the dimensionless parameters ^ 
and u = V jT as well as on the particle volume fraction. 
Unless otherwise stated, our simulations are done at vol- 
ume fraction 4% (i.e., ■kNg'^IIoV = 0.04), with iV^lOOO 
and ^=0.11(7, and we take oq — O.lScr. We quote times 
in MC steps (units of tq), noting that tb ~ 44ro for our 
chosen step size qq. We observe that this step size oq is 
comparable to the interaction range ^cr, so that our re- 
sults are not yet representative of the limit of small uq. 
We have conducted simulations with smaller oq: while 
quantitative differences are observed, qualitative features 
are unchanged. As usual, our time step (or equivalently, 
flo) is chosen as a trade-off between accuracy of numerical 
integration and practical efficiency. 



where rii is the number of neighbours for particle i: that 
is, the number of particles within the interaction range 
^cr. It is also convenient to define Ei — — ^riiU. 

This model is a simple but effective description for col- 
loids and some globular proteins [55J [3T1 ■ As shown 



III. CRYSTALLISATION PROCESS 

In Fig. [2] we show results from dynamical simulations 
for various bond strengths u. In the stable fluid phase 
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FIG. 2. (a,b) Average number of neighbours n{t), plotted as 
a function of time, for different bond strengtfis. Each curve is 
averaged over eight trajectories, (c) The 'crystaUinity' of the 
system as a function of time for the same six bond strengths. 
The definition of the order parameter ni42{t) is discussed in 
the text. Each curve is an average over eight trajectories. 

(e.g., u = 1.85), the system quickly equilibrates in a fluid 
state where the average number of neighbours per parti- 
cle 

n{t) = (n.(t)) (3) 

is relatively small. Increasing the bond strength to 
u — 2.3, the system is in the fluid-solid phase coexistence 
region of the phase diagram, and nucleation is observed 
at a time i ss 6 x 10^ MC steps. For stronger bonds, 



u = 2.6 and u = 3.3, the behaviour of n{t) shows that 
the nucleation barrier is small so that clusters of particles 
grow smoothly, starting at early times. For very strong 
bonds u > 4.5, clusters of particles grow rapidly at early 
times, but this growth slows down at longer times. This 
is the kinetic trapping regime. 

Extrapolating from the results of Liu et al. f28^, we 
believe that the state point u = 1.85 is in the dilute fluid 
phase (outside the fluid-crystal binodal) while the point 
u = 2.3 is close to the metastable liquid-vapour binodal. 
For simulations at u = 2.2 (not shown), no nucleation 
was observed for times up to 5 x 10^ MC steps. This is 
consistent with earlier observations pOl [55] that nucle- 
ation is rapid in the vicinity of the metastable binodal 
but slow in the region between crystal-fluid and fluid- 
fluid binodals. 

We use a common neighbour analysis (CNA) [29] to 
measure the crystaUinity of the assembling system. In 
particular, we count the number of bonded pairs of par- 
ticles that obey a 'crystaUinity criterion'. The criterion 
is that pairs of particles have exactly four mutual neigh- 
bours, and those mutual neighbours must share exactly 
two bonds. We denote the number of bonded pairs that 
satisfy this criterion by N142, and we normalise it by 
defining ni42(t) = 2{Nn2{t)) /N . This normalisation 
facilitates comparison with the number of neighbours 
n{t) plotted in Figs. [SJa) and (b): if all bonded pairs 
in the system satisfy our 'crystaUinity criterion' then 
^142(0 = n(t). In the notation of Honeycutt and An- 
dersen [29j , we are measuring the combined number of 
1421 and 1422 environments, which are indicative of cu- 
bic crystal structures. 

In Fig. [2][c) we show 71.142 {t) for the same six indicative 
bond strengths considered so far. At early times 71142 
is small, with a sudden increase at later times, as crys- 
tallites form in the system. Taking the data from this 
figure aXt = 10^ MC sweeps, we obtain the 'yield' shown 
in Fig. [1] (Plotting the yield at later times shows similar 
results and is discussed in Section [TV C[ ). 

To summarise the results of Fig. [2] for m < 1.85, the 
system does not assemble. We identify a range of good 
assembly between u=2.3 to 3.3 with the optimum around 
M=2.4 to 2.8. For energies m=4.5 and above we consider 
the system to be dominated by kinetic trapping. 

In Fig. [3]Ja), the number of bonds is plotted as a func- 
tion of time for various temperatures, including the ones 
we have been focussing on up to now- the colours and 
symbols of the data are now mixed. The graph has been 
cropped so as to show the evolution of the system only up 
to 4000 time steps. We pose the following question: if we 
are not prepared to wait for the late-time information dis- 
cussed earlier in this section, can we determine the fate of 
the system by looking at the dynamics so early on? More 
specifically, which systems will remain fluid, which ones 
will assemble well, and which ones will get trapped? In 
the remainder of the paper we demonstrate that the only 
information needed is, in fact, in the narrow window in- 
dicated by the dashed lines between 2000-4000 MC steps. 
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FIG. 3. (a) We plot n(t) for various bond strengths, includ- 
ing the ones shown in Fig. |2] with the colours and symbols 
mixed up, and showing only t < 4000 MC sweeps. From 
this information alone, we aim to predict which systems will 
remain fluid, which will assemble into crystalline structures 
and which ones will get trapped. The dashed lines define the 
narrow window over which the dynamical measurements will 
be made, (b) The simulation protocol used to measure the 
correlation and response functions, as described in the main 
text. 

We find that the long-term fate of the system is strongly 
correlated with certain early-time measurements that we 
will describe. 



IV. PREDICTING ASSEMBLY QUALITY 

A. Reversible bonding, correlation and response 
functions 

We aim to predict assembly quality, assuming that 
all system information is available for 2000 < t < 4000 
MC sweeps. The hypothesis is that the dynamics of the 
system at early times contains enough information on 
whether the system will assemble and how well. (The 
specific time range 2000 — 4000 MC sweeps is chosen 
only for concreteness: the dependence of our results on 
the choice of time window is discussed below.) A cen- 
tral task is to identify the propensity of the system for 
assembly or kinetic trapping. Kinetic traps arise through 
bonds which do not allow for long range crystalline order. 
Whitesides [35] observed that the reversible formation 
of weak bonds allows the system to escape from kinetic 
traps, facilitating good assembly. Conversely, if attrac- 
tive interactions between particles are too strong, then 



particles tend to aggregate into disordered clusters that 
grow rapidly and do not anneal into crystalline struc- 
tures. Our goal is to make measurements that exploit the 
relationship between bond-breaking and effective assem- 
bly, in order to distinguish which systems will assemble 
and which will suffer from kinetic trapping. 

To this end, we consider correlation and response func- 
tions that depend on the behaviour of the system be- 
tween two times. Away from equilibrium, correlation 
and response measurements have been investigated in 
the theoretical study of glassy systems [36 3^, in the- 
ories of non-equilibrium processes |39l HP] as well as in 
self- assembling [30j and gelating systems [41]. Moreover, 
experimental developments have shown that measuring 
correlation-response functions may be a possible route 
for characterising dynamics in the lab too [42H44] . 

We define the (dimensionless) single-particle energy 
autocorrelation function as 

Cit,t^) = l^[{E,{t^)E,{t)) - {E,{t^)){E,{t))] 

= \[{n,{t^)n,{t)) - {n,{t^)){n,{t))] (4) 

We adopt the convention that < 't■^ consistent with 
Fig. [3]jb). For fixed ^w, the correlation function mea- 
sures the extent to which the system's structure at t„ is 
correlated with its structure at some later time t. For 
example, in the equilibrium fluid state the C(t,t^) de- 
cays to zero on a time scale that reflects the lifetime of 
an interparticle bond. 

Away from equilibrium, it is convenient to normalise 
this correlation function as 

The equal time correlation function C{t,t) measures the 
variance in the number of bonds between particles, at 
time t. If bonding is irreversible (bonds never break), 
C{t,t^) may be interpreted as the fraction of bonds in 
the system at time t that had already been formed at 
the earlier time t^. Increasing t at constant t^, the sys- 
tem forms new bonds, and the correlation C(t, t^) decays 
towards zero. Thus, in contrast to the situation at equi- 
librium, the correlation function away from equilibrium 
does not simply measure a bond lifetime, but a combina- 
tion of a bond lifetime and the rate of bond formation. 

In order to separate these effects, we also consider a 
response function. The idea is that bond-making occurs 
whenever particles collide, so their rates are largely inde- 
pendent of the bond strength. On the other hand, bond- 
breaking processes require thermal activation over an en- 
ergy barrier U and take place with rates proportional to 
Q-u/ksT ^ On changing the strength of the forces between 
particles, the bond-breaking rate will change: measuring 
this response gives information about the relative likeli- 
hood of bond-making and bond-breaking in the system. 
The protocol for measuring this response is shown in 
Figl3l;b). As before, we begin with an equilibrated system 
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FIG. 4. Correlation C{t,t^) (solid lines) and response func- 
tions t-w) (symbols) are plotted versus time t for three dif- 
ferent bond strengths. We also plot 1 — C{t, t^) (dotted line), 
which would be equal to the response if the system were at 
equilibrium, (a) Bond strength ii=1.85: the system is equili- 
brated in a dilute fluid phase. The response is shown in red 
circles it follows the dotted line (i.e., x(i, iw) = 1 — C{t,t„) 
in accordance with FDT). (b,c) Bond strengths u = 2.6 and 
u = 4.5 respectively: both systems are far from equilibrium 
and the responses (symbols) differ from the plots of 1—C{t, t„) 
(dotted lines). For the weaker bonds (panel b), the system 
will assemble into a crystal; for the stronger bonds (panel c), 
the system will become kinetically trapped. 



of hard spheres, and introduce the interaction potential 
U at time t = 0. After a waiting time tw we perturb 
the bond strengths in the system, so that the energy of 
the ith particle becomes Ei = —^{U + 5Ui)ni, where 
Hi is the number of neighbours of particle i (as above), 
while 6Ui is the perturbation applied to the ith particle, 
and the factor of ^ ensures no double counting of bonds. 
At a later time t > ty,, we measure the (dimensionless) 
response of particle i to the change in its bond strength: 



d{n,{t)) T 



(6) 



This response depends on since the perturbation is 
applied only between times t and t^. In our simulations 
t and <w vary within the window of 2000-4000 time steps 
as indicated by the dashed lines in Fig.jsjja). 

As shown in Appendix |Aj x(t,iw) can be measured 
computationally by applying perturbations of magnitude 
5Ui to all particles, but with randomly chosen signs. 
The response is then obtained by comparing those par- 
ticles with positive and negative values of 5Ui. We 
take \5Ui\ = 0.125C/ and have checked that for all bond 
strengths considered the perturbation is within the linear 
response regime. 

In Appendix |Bj we prove the fluctuation-dissipation 
theorem (FDT) that links C{t,t^) and x(^:^w)- Normal- 
ising the response as 



C{t,t) 



the FDT reads 



(7) 



(8) 



where we added the label 'eqm' to emphasise that this 
relation holds only at equilibrium. 

Fig. ^a.) shows correlation and response functions that 
are typical for a system in the dilute fluid phase at equi- 
librium (the bond strength is m = 1.85). The perturba- 
tion is applied at tw — 2000 MC steps. The correlation 
function, is maximal at t = and then decays to zero, 
showing that the system eventually loses all memory of 
the bonds it made at time The response to the ap- 
plied perturbation, starts to grow immediately after the 
perturbation is applied. We also plot 1 — C{t,tyj) which 
is equal to the measured response, in accordance with 
FDT, since the system is at equilibrium. 

Assembling systems, however, are far from equilib- 
rium. Nevertheless, FDT can provide a useful compari- 
son quantifying how far the systems deviate from "locally 
equilibrated" states From the previous section we 

know that the system with u = 2.6 will assemble into a 
crystal, within ~ 10^ MC steps. In Fig. |4|^b) we show 
correlation and response functions at times of order 4000 
MC steps, long before crystal formation. The decay of 
the correlation is slower than for u = 1.85 and it does 
not reach zero within the range of times shown, indicat- 
ing that the system retains some memory of the bonds it 
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made at t^. It is also clear that the two curves x(*7^w) 
and 1-C'(i,tw) no longer trace each other. In Fig. jij^c) 
we show the same plot for u ~ 4.5. The correlation func- 
tion decays very slowly indeed. The response is almost 
zero and there is a large deviation between the response 
x{t, tw) and 1-C'(i, t„). A near-zero response implies min- 
imum unbinding and therefore little annealing of disor- 
dered bonds. The bonds that were made at iw persist for 
long times, throughout the 2000 time-step window shown 
here. This indicates that the system will eventually get 
kinetically trapped, as was indeed shown in the previous 
section. 

The significance of the above observations is that 
we have identified a measurement which distinguishes 
between different far-from-equilibrium regimes of be- 
haviour, using only information from much earlier times, 
at which the system bears little structural resemblance 
to its final product. Recalling Fig. [2jc), the measure of 
crystallinity 77,142 is tiny on these time scales {t ~ 10'^) for 
all the bond strengths considered. At these early times, 
a structural analysis of the system has little predictive 
power. On the other hand, the data of Fig. |4] change 
considerably as we pass between the equilibrated fluid 
phase (weak bonds), the good-assembly regime (interme- 
diate bond strength) and the kinetically trapped (strong 
bond) regime. 

B. Estimating fluctuation-dissipation ratios 

In the context of glassy model systems, the clearest 
way to analyse correlation-response data is to fix the time 
t and vary HD ESI ES] . Note that this is not the case 
in Fig. |4] where correlation and response data are plotted 
as a function of t at fixed tyj. Recalling Fig.[3](b), we now 
fix the time t at which the measurement is made and vary 
the time t„ at which the perturbation is switched on. 

Results for correlation and response functions are 
shown in Fig. [sja) and (b) respectively. For t — t^, 
the response vanishes since the perturbation has had no 
time to act; as decreases towards zero, there is an 
increase in the time t — over which the perturbation 
acts, so the response grows (going from right to left in 
Fig. [Sjjb)). The gradient dx/dtw has an interpretation 
as the response of the system to an instantaneous (im- 
pulse) perturbation applied at t^. The only difficulty 
when obtaining data with fixed t and variable is that 
each value of requires a separate computer simulation 
because the time at which the perturbation is applied is 
different in each case. 

In Fig. |6j we summarise the data of Fig. [5] by mak- 
ing a parametric plot of the response x(t,iw) as a func- 
tion of the correlation C{t,t^), keeping fixed t — 4000 
MC steps and varying and the bond strength. Such 
fluctuation-dissipation (FD) plots are often used in the 
study of glassy systems [36l IHl EH ESI to summarise 
measurements of correlation and response functions. The 
bottom right corner, where the correlation is maximum 
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FIG. 5. (a) Correlation C{t,t„) plotted against iw for dif- 
ferent bond strengths at fixed t = 4000. (The symbols for 
the characteristic bond strengths plotted here were defined in 
Fig. [ija) and are also shown in the legend of Fig. |6|. When 
t = the correlation is maximum and decays as decreases 
(or t — tw increases), (b) Response x(^)^w) plotted against 
for different bond strengths at t = 4000 MC steps. Correla- 
tion and response functions are estimated by averaging over 
many independent trajectories, with error bars shown for a 
few representative points in each panel. 



and the response zero, corresponds to t = tw Follow- 
ing the curves from right to left (decreasing C'(i,t„)), 
the points indicate the behaviour as decreases. The 
dashed line is Xoqm{t,tw) = 1 — Ccqmit^t^), which is the 
FDT prediction for a system at equilibrium. The data 
for the high temperature system lie on the FDT line, 
as expected since the system is equilibrated. On crossing 
the binodal, the bond strength is sufficient to drive phase 
separation and assembly. In this regime, the data lie on 
the FDT line when t — t„ is small, before deviating when 
C and iw get smaller (see for example u — 2.3,2.6,3.3). 
The time window over which the data lie close to the 
FDT line decreases as the bonds become stronger. The 
low temperature data, u = 4.5, 7.0, deviate from the FDT 
line even when t — is very small. 
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FIG. 6. Parametric plot of x(i,iw) versus C{t,t„). The data 
are plotted for fixed t=AQOQ MC steps, varying between 
2000-4000 MC steps, in steps of 100. The dashed line is 
x(*i ^uj)=l"C'(t, tu,), which is the FDT behaviour at equilib- 
rium. Deviation of the data from the FDT line indicate devi- 
ations from local equilibration (see main text). 



Jack et al. [30] , argued that if the data are close to the 
FDT line when t — k, t , then the behaviour of the 
system is 'locally equilibrated' Wf\ on the time scale r. 
In this case, the idea is that the region of configuration 
space explored on these time scales is being explored re- 
versibly. That is, given a movie of the system of length 
T and a similar movie where time has been reversed, it 
would be difficult to discern which is which. Clearly, if 
this property holds, then the likelihoods of bond-making 
and bond-breaking must be similar. On long time scales 
(large t — t^), it will be apparent if a system is assem- 
bling and bond-making dominates. But if the system is 
avoiding kinetic traps by exploiting the reversibility of 
the bonding process in order to anneal out defects ^35j 
then one expects the system to be locally equilibrated 
over a time interval r that is comparable to the bond 
lifetime, and hence that the data in Fig. [6] remain close 
to the FDT line over a significant range of C. As in 
Ref. 1^ , this expectation seems to be borne out by these 
numerical simulations. 



C. Robustness of results 

Comparing Figs. [T] and |6j long-time measurements of 
the yield are correlated with short-time measurements 
of correlation and response functions. We have in mind 
that the short-time measurements might be used to pre- 
dict the long-time behaviour. However, these systems are 
far-from-equilibrium, and both the short- and long-time 
measurements will depend on the times at which these 
measurements are made. If the short-time measurements 
are to be a useful predictive tool, the correlation between 
short-time and long-time behaviour must be robust to 




-- FDTline 
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0.4 0.6 

C(t,tJ 



FIG. 7. (a) Yield at various times, plotted as a function of 
the bond strength (compare Fig. [TJa)). (b) Parametric plot 
of x(^j^w) versus C{t,t„) for different measurement times: 
t = 4000 MC steps and 2000 < < 4000 MC steps as in 
Fig. [6] and t = 12000 MC steps with 6000 < iw < 12000 MC 
steps. The dashed line is x{^y'tw)=^-C{t,tw), which is the 
FDT result at equilibrium. 



variations in the time at which the measurements are 
made, as well as to changes in the system parameters 
(for example, volume fraction and interaction range). 

In Fig. [T) we show the effect of varying the mea- 
surement times, keeping system parameters constant. 
Fig. [7](a) shows that measuring the yield at different 
times leads to differences in the crystallinity of the sam- 
ple, as expected since the phase transformation is taking 
place over the whole time window considered. Neverthe- 
less, a change of two orders of magnitude in the mea- 
surement time leads to the same qualitative results, and 
the condition that assembly is optimal for u ss 2.5 is 
robust. In Fig. ^h), we show that increasing the time 
at which the correlation and response measurements are 
made leads to very small changes in the FD plot. (We 
show results for the near-optimal condition m = 2.6 but 
results are similar for other bond strengths. Compared 
to Fig. [6j we have increased t and by a factor of 3.) 
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Changing ^ from 0.11 to 0.25 has a significant effect on 
the yield (crystallinity) of the self-assembly process. In 
particular, for ^ = 0.25 and a fixed time t — 10® MC 
sweeps, Fig. [sj^c) shows that the yield rather is small for 
all u. We find that on increasing the measurement time, 
the yield increases significantly, but the position of max- 
imal yield remains stable at u 2 (data not shown). In 
Fig.jsj^d), we show FD plots for various ^, where the value 
of u is chosen in each case to be near near-optimal for the 
long-time yield. Quantitative differences are apparent as 
^ is varied, but we believe that the FD plots at optimal 
assembly are similar enough for predictions of effective 
crystallisation conditions to be made. 



OUTLOOK 



FIG. 8. (a) Measurement of the yield at t = lO'' MC steps, 
as a function of the bond strength, for three different vol- 
ume fractions, (b) FD plot, calculated at optimal assembly 
u = 2.6, for two three different volume fractions, (c) Yield 
measurements for three different interaction ranges, (d) FD 
plot, calculated using the optimal value of u associated with 
each value of ^. While the value of it changes, the FD plot for 
optimal It is qualitatively unchanged. 



This insensitivity to changes in t and tw reinforces the 
idea that the FD plot has potential as a predictive tool, 
since one arrives at the same prediction, regardless of the 
specific time at which the measurements are made. 

We now discuss dependence on the volume fraction 
and the interaction range ^. Results are summarised in 
Fig. [8| On increasing the volume fraction from 4% to 
16%, the yield (crystallinity at long times) varies only 
slightly. Increasing at constant u increases the super- 
saturation and hence speeds up nucleation, but this ef- 
fect is rather weak in this system, for these times. In the 
simulations of Fortini et al. |26j . nucleation was found 
to be very slow for values of u between the crystal-fluid 
and fluid-fluid binodals, and they observed crystallisa- 
tion only when metastable fluid-fluid phase coexistence 
was also possible. Our results are consistent with this ob- 
servation, although we were not generally able to iden- 
tify long-lived metastable phase coexistence in our dy- 
namical simulations. In the kinetic-trapping regime, the 
dependence of the results on volume fraction are surpris- 
ingly weak: this is presumably related to the thermally- 
activated time scales for bond-breaking, which are inde- 
pendent of (j). Turning to the correlation-response mea- 
surements in|8][b), the dependence on volume fraction is 
again weak, providing further evidence that FD plots can 
be used to predict conditions for effective crystallisation. 

On increasing the range of the attractive interactions, 
the phase diagram of the system changes as the liquid 
phase becomes more stable [28!. For the case ^ = 0.25, 
the fluid-fluid and crystal-fluid binodals are very close 
to each other (both are near to w = 1.5 for = 4%). 



We have considered how we might predict the evolu- 
tion of a model colloidal system based on the dynamics at 
early times, as presented schematically in Fig. [3](a). The 
data shown in Figs. [6] [7] and |8] seem promising in this re- 
gard. Such parametric plots can be used during the early 
part of the crystallisation process to distinguish between 
different assembly regimes, that will become structurally 
apparent only at later times. The current paper comple- 
ments a previous study 30J, indicating that the method 
is applicable at least to patchy particles that form viral 
capsids, as well as to discs and spheres with short-ranged 
attractive interactions. These systems have diverse kinet- 
ically trapped states and may have no equilibrium phase 
transitions at all (viral capsids), or phase transitions to 
both stable crystalline and metastable liquid states (as 
in this work). Despite these differences, the correlations 
between parametric plots such as Fig. [6] and yield mea- 
surements such as Fig.[IJb) seem to be conserved between 
models, indicating that the method may have broad ap- 
plication. 

We believe that the qualitative similarity between sys- 
tems is due in part to the observation of Whitesides [3S] 
that reversibility is essential for effective assembly, in a 
wide variety of contexts. Essentially, for a system to 
assemble well into an ordered structure it has to have 
the ability to undo incorrect bonds, and thus avoid ki- 
netic traps [30l |35l HH |49j . Excessively strong interac- 
tions lead to irreversible sticking and produce disordered 
states: this is avoided only if assembling components can 
adjust their positions even after they have already formed 
bonds. Whitesides |35] referred to this property as re- 
versibility and stated it as a qualitative requirement for 
good self-assembly or crystallisation. FD plots such as 
Fig. 6 constitute one route towards a quantitative ver- 
ification of this statement. When t — t^ is small (and 
C is large), the system is close to local equilibrium as 
discussed in Section |IVB[ the range oit — t^ over which 
local equilibration holds is an interval of reversibility r. 
Fig. [5] shows that for systems that assemble effectively 
then r is relatively large, while r tends to be small for 
systems that eventually get arrested. 
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Further, the parametric plot in Fig. [6] provides a 'di- 
mensionless' criterion for separating these regimes: in- 
stead of the time r over which the system behaves re- 
versibly, it is natural to identify a range of C, which is 
a dimensionless variable that can be compared directly 
between different systems (for example, one may com- 
pare the results of this article with those of Ref. '30'). 
The only free parameter associated with this compari- 
son is the time t associated with the parametric plot: we 
have also shown in Fig. [T^a) that results for this system 
depend weakly on this time, strengthening the argument 
that comparison of plots from different systems under dif- 
ferent conditions can be compared fairly with each other. 
In the future, we hope that the predictive information in 
such measurements might be useful in optimising com- 
puter simulations of assembling systems, and perhaps 
even experimental crystallisation and self-assembly pro- 
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Appendix A: Calculation of response function 

This section explains how the response function in 
Eq.(|6| is calculated in simulations. For compactness of 
notation, we write hi — 5Ui/2T so the response of Eq.([6| 
is 



4 dh, 



(A.l) 



where the field is applied between times t and i„. We 
run simulations where hi is finite for all particles: we 
take hi — \h\ for half of the particles (chosen at random), 
with hi — — \h\ for the other half. 

In the presence of such a perturbation, we make a Tay- 
lor expansion of {ni)h, where the subscript h indicates 
the presence of the N perturbing fields hi, ... , h^. The 
result is 



therefore write, for i ^ j. 



djui 
dhi 



c 

N 



(A.3) 



with c = 0(1) as iV — oo (see below). 



We now decompose the sum over j in (A. 2) into a 



contribution from those particles j for which hj > 0, and 
those with hj < 0. Restricting to j ^ i, let 5*+ be the 
set of particles j for which hj > and S- the set with 



hj < 0. In (A.2), all terms in the sum are equal to either 



-r\h\c/N or ~\h\c/N, so if the number of particles in S+ 
is 7V+ and the number in S- is N- then we have 



djui) 



\h\c 
N ■ 



(A.4) 



and hence 



{ni)h = (rii) + h, 



djui 
dh. 



\h[ 



' N 



0{h^) (A.5) 



It is clear from this equation that we require c = 0(1) as 
N oo (see above) so that the response is finite if (for 
example) a positive field is applied to exactly half of the 
particles and the other particles are unchanged = ^ 
and iV_ = 0). 

In the case where half of the particles receive a pertur- 
bation of +\h\ and half receive — the value of {ni)^ 
depends on i only through the sign of hi. For particles 
with hi > 0, we write their average number of neighbours 
by {ni)j^ while for particles with /i^ < we write {rti)^. 
These quantities are readily calculated in a simulation, 
by evaluating the number of bonds for each particle and 
averaging separately over those with hi > Q and those 
with hi < 0. 

In total, ^ particles have hi > and ^ have hi < 0. 
However, the sets 5+ and S- both exclude particle i, so 
if /ij > then 7V+ = f - 1 while iV- = f . Hence 



\h\ 



djui 
dh. 



\h\ 



N 



0{h^ 



(A.6) 



Similarly if h, < 0, then N+ = ^ while = f - 1, so 
that 



\h\ 



dh. 



\h\- + 0{h') 



(A.7) 



{ni)h = {n-i) + h, 



dhi 



^ dh, 



+ Oih^) (A.2) 



where all derivatives are evaluated ai h — and we omit 
the dependence of rii on time t, for brevity. The first- 
order terms in the Taylor expansion are the response of 
particle i to its own field hi, and the response of this 
particle to the specific combination of other fields. 

For j 7^ i, the response is independent of i and 

j, and scales as in the thermodynamic limit. We 



It is therefore clear that the average particle energy 
in the absence of the field may be estimated in the per- 
turbed system by 

{n.)^^^^^^^±±^+Oih^) (A.8) 

while the single particle response function can be esti- 
mated as 



djui 
dh. 



2\h\ 



+ 0i\h\) + 0il/N) (A.9) 
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That is, the response ^w) in Eqs.(l6|) and ( A.l I can be 
estimated as iw) = 4\sij\ i{n'i) + ~v^i)-) by measuring 
the difference in the number of neighbours of particles for 
which the perturbing field SUi has positive or negative 
sign. This is the method used to calculate x(ij ^w) hi this 
article. 



In (B.l), the only /i-dependence comes through G'' so 
we seek an expression for g^G'' . We use detailed balance 



together with (B.3) to write 



Gl,(I^F) = Gl,(F^iy. 



^uB(F)-uB(I)+Y,i[hini(F)-hin,(I)] 



(B.6) 



Appendix B: Fluctuation-dissipation theorem 

In this section, we prove Eq.Q of the main text. As 
in the previous section, we write hi — bUijlT. The re- 
sponse of Eq.([6| may be written in terms of probabilities 
by using (A.l) together with 



d_ 
dhi 



t-tjl 



F)n,{F) , (B.l) 



I.F 



where / is the configuration of the system at time and 
F is the configuration at time t: the sum runs over all 
possible configurations / and F. Also, ni{F) is the num- 
ber of neighbours of particle i in configuration F] the ini- 
tial distribution p{I) is the probability of being in config- 
uration / at time t^,; and the propagator Gf_f^{I — ^ F) 
is the probability of being in configuration F at time t, 
given that the system was in configuration / at time t^. 
The label h on the propagator indicates that it depends 
on the applied fields hi, while p{I) and ni{F) do not. 

Since we are concerned with Eq.Q, we restrict to the 
case of equilibrium response functions, for which the sys- 
tem is equilibrated with /i — at time t^: 



p{i) = pt^Ai) = 



^uB(I) 

z 



(B.2) 



where ^ = E/ e"^^^^ and B{I) = \ Ml) is the total 
number of bonds in configuration /. The energy of con- 
figuration / in the unperturbed system is —UB{I) while 
in the perturbed system it is —UB{I) — T^^hini{I). 
Hence we define 



h 

r cqm 



{I) 



1 



,nS(/)+5:,hin,(/) 



(B.3) 



with Zh = Y^jc'^BW+T..h,n,ii) ^j^g partition function. 

To prove Eq.([8| we make use of detailed balance. For 
a single Monte Carlo step in the presence of the pertur- 
bation h, let the probability of arriving in configuration 
J from an initial configuration / be P^{I — )■ J). Detailed 
balance states that 



^cqm 



{i)p\i ^ J) = p^^,M)pH J ^ I) 



(B.4) 



which ensures that the system converges to the equilib- 
rium distribution in the limit of long times. A similar 
relation follows for the propagator: 

pLmiI)Gt-t„ il^F)^ Pe\m(^)Gt. (F ^ /). (B.5) 



Taking a derivative with respect to hi and evaluating it 
at h — 0, we arrive at 



dh. 



(n,(F) - umG^tltiF + Qj^Gl,jF ^ /) 



Now, 



detailed 



balance implies 



(B.7) 



that 



^uB{F)-uB(i)Qh^o^l^p -> /) = G'lziil ^ F), and 
we also have pfj=/(i^)e^)~"-^(^2^ p"cq=m(^)- Combin- 



ing these results with (B.l) and (B.7 1 yields 
d{ni 



dh^ 



E-^otmWGjT^J/ ^ F)[n,{F) - n.(/)]n.(F) 



d 



I,F 



where we recognise the first term on the right hand side as 
the correlation function {ni{t)[ni{t) — ni(iw)])7 evaluated 
at equilibrium. The second term on the right hand side 
is zero since '-^J'-tw ^ ^) = 1: this follows from the 
definition of Gt-t^ {F — ^ /) as the probability of being in 
state I at time t since these probabilities must sum to 
unity, regardless oi F, h and t — t^. 



Hence, at equilibrium, (B.8) reduces to 

d{ni) 



dh. 



{ni{t)[ni{t) - n^(tw)])- 



(B.9) 



To recover Eq.([8| of the main text, we note from ^ 
that 

G(i,iw) = \[{n,{t)n,{t^)) - {n,{t)){n,{t^))] (B.IO) 



so that the right hand side of (|R9]) is 4[G(t, t)-C{t, + 
{ni{t)){ni(t) — ni{tyj)). At equilibrium, time-translational 
invariance implies {uiit) — riiit^)) = 0, and using (A.l) 



with (B.9) one arrives at the fiuctuation-dissipation the- 
orem 



Xcqin(^7^w) ^cqin(^iO ^cqm(^j^w)- 



(B.ll) 



Finally, dividing both sides by Ccqm{t,t) yields the nor- 
malised version of the FDT, which is given in Eq.(l8|). 
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